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Abstract 

In this paper, we investigate the emergence of a predator-prey system with Ivlev- 
type functional response and reaction-diffusion. We study how diffusion affects the 
stability of predator-prey coexistence equilibrium and derive the conditions for Hopf 
and Turing bifurcation in the spatial domain. Based on the bifurcation analysis, we 
give the spatial pattern formation, the evolution process of the system near the 
coexistence equilibrium point, via numerical simulation. We find that pure Hopf 
instability leads to the formation of spiral patterns and pure Turing instability 
destroys the spiral pattern and leads to the formation of chaotic spatial pattern. 
Furthermore, we perform three categories of initial perturbations which predators 
are introduced in a small domain to the coexistence equilibrium point to illustrate 
the emergence of spatiotemporal patterns, we also find that in the beginning of 
evolution of the spatial pattern, the special initial conditions have an effect on the 
formation of spatial patterns, though the effect is less and less with the more and 
more iterations. This indicates that for prey-dependent type predator-prey model, 
pattern formations do depend on the initial conditions, while for predator-dependent 
type they do not. Our results show that modeling by reaction-diffusion equations 
is an appropriate tool for investigating fundamental mechanisms of complex spa- 
tiotemporal dynamics. 
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1 Introduction 



A fundamental goal of theoretical ecology is to understand how the inter- 
actions of individual organisms with each other and with the environment 
determine the distribution of populations and the structure of communities. 
Empirical evidence suggests that the s patial scale and stru c ture o f environment 
can influence population interactions (ICantrell &: Cosnerl . |2003|). The endless 
array of patterns and shapes in nature has long been a source of joy and wonder 
to laymen and scientists alike. Discovering how such patterns emerge sponta- 
neously from an orderless and homogeneous environment has been a challenge 
to res earchers in the natural sciences throughout the ages (IBen- Jacob fc Levind . 
20011 ) . The problem of pattern and scale is the central problem in ecology, uni- 
fying populat ion biology a nd ecosystems science, and marrying basic and ap- 
plied ecology (ILevinl . Il992l ) . The study of spat i al patterns in the distribu t ion of 



organisms i s a central issue in ecoloev (iLevid. [19921; iKoch fc Meinhardtl , 
Neuhauserl.l200ll:lAlonso. Bartumeus &: Catalan. 2002: Medvinskv et al . 



Cantrell fc Cosnei . 2003 : Leppanen et al . 20031: Leppanen . 2004; Murrav 



Hawick. James &: Scogingsl.l2006l:lGrifEth &: Peres-Netobl.l2006l:lMaini. Baker ^ Chuong 




20061: IPearce et aj |2006|; iBaurmann. Gross, fc Feudel l2007l: iLiu fc Jinl. hom. 

Schnell. Grima fc Main! 200?! : Shoji. Yamada. Uevama fc Ohta . 200?! : Wang. Liu fc Jin 
20071 ) since the pioneering work of Alan Turing (jTuingl . Il952l ). And Turing 
reaction-diffusion system exp lains spatial pattern s spontaneously forming in a 
perfectly homogeneous field (lUriu &: Iwasal . 120071 ). The instability now identi- 
fied with Turing's name is believed to be involved in th e formation of structure 
in many systems of biological interest (iMurrayl . |2003| ). Theoretical work has 
shown that spatial and temporal pattern formation can play a very important 
role in ecological and evolutionary systems. Patterns can affect, for example, 
stability of ecosystems, the coexistence of species, invasion of mutants and 
chaos. Moreover, the patterns themselves may interact, leading to selection 
on the level of patterns, i nterlocking ecoevol utionary time scales, evolutionary 
stagnation and diversity (jSavill fc Hogewegl . [T999i l 



The or i gin of these patterns has commonly attributed to two sorts of sources 
( jLevinl . Il992l ). One is a heterogeneous distribution of abiotic factors and the 
other is underlying mechanisms at the level of individuals. Patterns generated 
in abiotically homogeneous environments are particularly interesting because 
they require an explanation based on the individual behavior of organisms. 
They are commonly called "emergent patterns", because they emerge from 
interactions in spatial scales that are much larger th an the characteristic scale 
of individuals (lAlonso. Bartumeus fc Catalanl . l2002l ). The instability leads to a 
process that might be called differentiation and in its simplest realization is the 
result of a competition between an activator (for predator-prey system, prey) 
and an inhibitor (for predator-prey system, predator) diffusing at different 
rates. The results of the instability has one characteristic property: its scale 
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or wavelength is determined by the concentrations of ambient species and the 
diffusion coefficients, and is therefore independent of any externally imposed 
length scales. In the process of morphogenesis the instability is likely to be 
triggered by the increasing scale of the system: the instability occurs once 
the system is large enough that it conta ins several natural wavelengths of the 
instability (jCallahan &: Knoblochl . Il999l ). 



The past investigations have revealed that spatial inhomogeneities like the in- 
homogeneous distribution of nutrients as well as interactions on spatial scales 
like migration c an have an important i mpact on the dynamics of ecologi- 



cal populations (iMedvinsky et all . l2002l : iMurrayl . |2003| ). In particular it has 



been shown that spatial inhomogeneities promote the persistence of ecologi- 
cal po pulations, play an important role in speciation and stabilize population 
levels ( jBaurmann. Gross, fc Feudell . 120071 ). Spatial ecology today is still domi- 
nated by theoretical investigations, and empirical studies that explore the role 
of space are becoming more common du e to technological advances that allow 
the recording of exact spatial locations (INeuhauserl . I2OOII ). 



On the other hand, as we know, our ecological environment is a huge and highly 
complex system. This complexity arises in part from the diversity of bi ologi- 



cal sp ecies, and also from the complexity of every individual organism (iJost 
19981 ). The relationship between predators and their prey has long been and 



will continue to be one of dominant themes in both ecolo gy and mathemati- 
cal e cology due to its universal existence and importance (IKuang fc Berettal. 



1998 



) ■ A classical predator-prey syste r n can be written as the form (lAbram fc Ginzburgl . 



2OOOI : lAlonso. Bartumeus fc Catalanl . l2002f ): 



it = uf{u) — vg{u,v), i) = h[g{u,v),v]v. 



:ii 



where u and v are prey and predator density, respectively, f{u) the prey 
growth rate, g{u,v) the functional response, the prey consumption rate by an 
average single predator, which obviously increases with the prey consumption 
rate, and can be influenced by the predator density, h[g{u, v),v] the per capita 
growth rate of predators (also known as the "predator numerical response"). 
The most widely accepted assumption for the numerical response is the linear 
one ( Arditi &: Ginzbum 19891 : Alonso. Bartumeus &: Catalan . 2002): 



h[g{u,v),v] = eg{u,v) - (3, 

where /3 is a per capita predator death rate and e the conversion efficiency of 
food into offsprings. 

In population dynamics, a functional response g{u,v) of the predator to the 
prey density refers to the change in the dens ity of prey attached p er unit time 
per predator as the prey density changes (IRuan fc Xiad . I2OOII ). In general, 
functional response can be classified as (i) prey dependent, when prey density 
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alone determines the response, g{u,v) = p{N); (ii) predator dependent, when 
both predator and prey populations affect the response. Particularly, when 
g{u,v) = p{^), we call model ([T]) strictly ratio dependent; and (iii) multi- 
species dependent, when species other t han the focal predator and its prey 
species influence the functional response ( lAbram fc Ginzburd . |2000| ). 



Th ere have been s everal famous function al response types: Hol ling type I- 
III flHollind . [l959al Jbl): Hassell-V arlev tvpe dHassell fc Varlevl .ll969h: Beddington- 
DeAngelis type by Beddington (IBeddingtonl. Il975[) and DeAngelis, Goldstein 
and Neill (IDeAngelis. Goldstein fc Neilll . Il975l ) independently; the Crowley- 
Martin tvpe fICrowlev fc Martini. Il989h: and the recent well-known ratio- dependence 
type by Ardit i and Ginzburg (lArditi &: Ginzburd . Il989l ) later studied by Kuang 
and Beretta ( Kuang fc Beretta . 19981 ). Of them, the HoUing type I — III was 
labeled "prey- dependent" and the other types that consider the interference 
among predators were labeled "predator-dependent" (lArditi fc Ginzburgj . ll989l ). 



Besides Holling type I-III, there is another impor tant prey-de pendent func- 



^ Hi 

tional response-Ivlev-type, originally due to Ivlev (llvlevl . Il961h : 



g{u,v) = l-e-^\ (2) 
and the corresponding Ivlev-type predator-prey model takes the form: 



u = u{l-u)-v{l-e~^''), v = ev{l-e-^'')-Pv. (3) 

where u and v represent population density of prey and predator at time 
t, respectively, e, P, 7 are positive constants, e the conversion rate of prey 
captured by predator, /? the deathrate of predator, and 7 the efficiency of 
predator capture of prey. From an ecological viewpoint, the conditions m > 
and V > must hold. From the second equation of ([3]), we know e > [3. 

In this paper, we mainly focus on the following spatial Ivlev-type predator- 
prey model with reaction diffusion: 

growth due to prey mortality due to prey random motility 

u= uiX^u) - 'v{l - e-^") + = f{u,v) + diV^u, 

growth due to predator mortality due to predator random motility 

= g{u,v) + d2Vh. 

where a = | > 1, and di, d2 are the diffusion coefficients of prey and predator, 
respectively, = + is the usual Laplacian operator in two-dimensional 
space. 
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Both ecologists and mathematicians are interested in the Ivlev-type pre dator- 
prey model and much progress has been seen in the study of model (l3 l) (lMa\ 



2007 



1981 



Metz fc Diekmannl. Il986l : iKooiil . 1 19961 : ISudd. Il998l: iTianl. l2006l: IWan 



Preedv et al. . 2006) and model (IH) (Sherratt. Lewis fc Fowlerl. 



Sherratt. Eagan &: L ewis. 1 9971: iKav &: Sherrattl . l2000l : iPearce et all . l2006l : lGarvie 



20071 : iPreedv et al 



19951: 



i2006i : lUriu &: Iwasal . 120071 ) . The results indicate that the 
-prey model (131) and (HI) have widely applicabilities in ecol- 



ogy, such as dynamics in predator-prey system 


(May. 


1981: Metz & Diekmann. 


19861: Sherratt. Lewis & Fowler. 19951: Sherratt. Eaean & Lewisl. 


19971: Tian. 


20061). host-Darasitoid system (Pearce et al. ! 


20061; 


Preedv et al.. 


2006), fish 


skin pattern (Uriu & Iwasa. 2007). and so on. 



Of them, Sherratt and co-workers ha d studied the dynamics of oscillations and 

chaos behind predator-prey inva sion (Sherratt. Lewis fc Fowleil. 1995 : Sherratt. Eagan fc Lewis . 
1997). Especially, in reference ( Sherratt. Eagan fc Lewisl . 1997 ). the authors 



performed a large number of numerical simulations of the invasion of prey 
by predator with four categories models, involving one- and two-dimensional 
reaction-diffusion model (jl]). For model (jl]), they used a large spatial domain, 
with the system initially in the prey-only steady state, except for a small region 
in the center of the domain, where a small density of predators was introduced. 
They stopped their simulations before the invading wave reached the end of 
the domain, so that the results were not sensitive to the boundary conditions, 
which could be either zero flux, periodic, or with population levels flxed at the 
prey-only steady state. They also discussed the way in which the populations 
evolved after the invasion had reached the edge of the domain. Furthermore, 
the authors performed a number of one-dimensional spatial patterns with the 
model (jl]) initially in the prey-only steady state and two-dimensional spatial 
patterns with two categories of initial perturbation. One is the introduction 
of predators are in a small, localized region of the domain, which is otherwise 
in the prey-only steady state; the other is the introduction of predators along 
a line running parallel to one edge of the (rectangular) spatial domain. Based 
on these results, the authors indicated that for model (j4]), the behavior be- 
hind the invasive front of predators consists of either irregular spatiotemporal 
oscillations, or periodic waves in population density. 



But to our knowledge, for model (jlj), the research on symbolic conditions 
of Hopf and Turing bifurcation, the evolution process of the spatial pattern 
formation, the mechanism of pattern formation emergence, especially the in- 
fluences of the speciflc choice of the initial conditions to the pattern formation, 
seems rare. 



The paper is organized as follows: In Section 2, we employ the method of 
stability analysis to derive the symbolic conditions for Hopf and Turing bifur- 
cation in the spatial domain. Based on these conditions we locate the Hopf 
and Turing bifurcation within the generalized parameter domain in 7 — c^i 
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bifurcation diagram. In Section 3, we give the spatial pattern formation, the 
evolution process of the system near the coexistence equilibrium point, via 
numerical simulation. For the sake of learning the influences of the initial 
conditions to pattern formation, we perform three categories of initial per- 
turbations which predators are introduced in a small, localized region of the 
circle, line and pitchfork domain to the coexistence equilibrium point to illus- 
trate the emergence of spatiotemporal patterns. Then, in the last section, we 
give some discussions and remarks. 



2 Stability and Bifurcation analysis 

The non-spatial model has at most three equilibria (stationary states), 
which correspond to spatially homogeneous equilibria of the model (jlj), in the 
positive quadrant: 

(i) (0, 0) (total extinct) is a saddle point; 

(ii) (1, 0) (extinct of the predator, or prey-only) is a saddle when 7 > — In 
or stable node when 7 < — In or saddle- node when 7 = — In 

(iii) a nontrivial stationary state {u*,v*) (coexistence of prey and predator), 
where 

= In(^) ' v* = -f In(^) (7 + ln(^)) = (5) 

with a > 1, 7 > — In 

In this paper, we mainly focus on the dynamics of nontrivial stationary state 
{u*,v*). For cyclical populations, this coexistence state will also be unstable 
and will lie inside a stable limit cycle in the kinetic phase plane. 

To perform a linear stability analysis, we linearize the dynamic system (jlj) 
around the spatially homogenous fixed point ([5]) for small space- and time- 
dependent fluctuations and expand them in Fourier space 

u{x,t) ~ u*e^^e^'''^, v{x,t) ~ ^*gAfgifc-x_ 

Then, in the linearized version of model (j4]), yielding a dispersion relation 
from which one can choose parameters to allow only some of the modes with 
Re(A) > to grow in time. The dispersion relation X{k) relating the temporal 
growth rate to the spatial wave number k can be found from the characteristic 
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polynomial of the original problem (jl]): 

A' - tr^A + Afc = 0, (6) 

where 

trfc = fu + gv- k'^idi + 62) = tiQ - k'^{di + ^2), 

Afc = fugv - fv9u - k^{fud2 + gvdi) + k'^didi = Aq - Pifnd2 + Qpdi) + k'^did2, 

here tro = fu + 9v, Aq = fuQv - fv9u, and the elements /„, fv,gu,gv are the 
partial derivatives of the reaction kinetics f{u,v) and g{u,v) denoted by (jl]), 
evaluated at the stationary state {u*,v*). 

The reaction-diffusion systems have led to the characterization of two ba- 
sic types of symmetry-breaking bifurcations — Hopf and Turing bifurcation, 
responsi ble for the emergence of spatiotemporal p atterns. See, for details, ref- 



erences flYang et all . 12002 : IWang. Liu fc Jinl . 120071 ). 



The Hopf bifurcation occurs when Im(A(A;)) 7^ and Re(A(A;)) = at A; = 0, 
and the critical value of Hopf bifurcation parameter 7 equals 

(1-«)(1 + A)' 



where A = (a — l)ln 



a-l 



At the Hopf bifurcation threshold, the temporal symmetry of the system is 
broken and gives rise to uniform oscillations in space and periodic oscillations 
in time with the frequency 

PA(-f{a - 1) + a) 
UH = lm{X{k)) = JAo = - ^ ^ 



7(a-l) 

and the corresponding wavelength is 

2n^^jia-l)^ (9) 



The Turing instability is dependent not upon the geometry of the system 
but only upon the reaction rates and diffusion. And it cannot be expected 
when the diffusion term is absent and it can occur only when the activator 
(for predator-prey system, prey) diffuses more slowly than the inhibitor (for 
predator-prey system, predator). Mathematically speaking, as di d2, the 
Turing bifurcation occurs when Im(A(A;)) = and Re(A(A;)) = at k = kT ^ 0, 



7 



kr is called the wavenumber. The critical value of bifurcation parameter 7 
equals 

_ A(l3A+2d2k^+d2Ak^) /-.qx 



where 



I Aq 
did2 



And at the Turing threshold, the spatial symmetry of the system is broken 
and the patterns are stationary in time and oscillatory in space with the 
wavelength 



A, 



27r _ 27r^rfici27^(a-l)2 
kr ~ J(3A{^{a-l) + A) 
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Fig. 1. (A) 7 — di Bifurcation diagram for the system (|4|) with a = 1.1, /3 = 0.5, 
(i2 = 0.2. Hopf and Turing bifurcation lines separate the parameter space into three 
domains. The other parameters in figures (B)-(F): (ii=0.02, and (B) 7 = 5.0; (C) 
7 = 5.552147102, the critical value of 7^; (D) 7 = 6.0; (E) 7 = 7.265163898, the 
critical value of 7^; (F) 7 = 8.0. The real and the imaginary parts of A(A;) are shown 
by solid lines and dotted lines, respectively. 

Linear stability analysis of model (jl]) yields the bifurcation diagram, the re- 
lation between 7 and di, is shown in figure [D^A). In this case, the Hopf bi- 
furcation line and the Turing bifurcation line separate the parametric space 
into three distinct domains. In domain I, located below all Hopf and Turing 
bifurcation lines, the steady state is the only stable solution of the system. 
Domain II is the region of pure Hopf instability, and in domain HI, located 
above all two bifurcation lines, both Hopf and Turing instability occur. 
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From the definition of Hopf and Turing bifurcation, we know tliat the relation 
between the real, the imaginary parts of the eigenvalue X{k) determine the 
bifurcation type. The relation between Re(A(A;)), Im(A(A;)) and k are shown in 
figure [T]^B)-(F). Figure [D^B) illustrate the case of parameter locate in domain I 
in figuredl^A), 7 = 5.0, one can see that Re(A(fc)) < and Im(A(fc)) 7^ at /c = 
0. Figure [It^C), 7 = 5.552147102 = •jh, the critical value of Hopf bifurcation, 
in this case, Re(A(fc)) = at = while Im(A(/c)) 7^ 0. In figure [11(D), 
7 = 6.0, the parameter locate in domain II, the pure Hopf instability occurs, 
one can see that at = 0, Re(A(A;)) > 0, Im(A(A;)) 7^ 0. Figure [11(E), 7 = 
7.265163898 = ■jt, the critical value of Turing bifurcation, at k = kx = 
2.116874108, Re(A(fc)) = Im(A(fc)) = 0. When 7 = 8.0, parameter locate in 
domain HI, figure [H^F) indicate that at /c = 0, Re(A(/c)) > 0, Im(A(fc)) ^ 0. 



3 Pattern formation analysis 



In this section, we have performed extensive numerical simulations of the spa- 
tially extended model in two-dimensional space, and the qualitative results 
are shown here. Model (jl]) is posed on a given domain Q = 400 x 400, with 
smooth boundary dQ. Zero-fiux Neumann boundary conditions are imposed 
on dfl to close the system. Model (jlj) is solved numerically with a = 1.1, 
P = 0.5, di = 0.02, d2 = 0.2 in two-dimensional space using a finite difference 
approximation fo r the spatial d erivatives and an explicit Euler method for the 



time integration (iGarviel . 120071 ) with a time stepsize of At = 1/3 and space 
stepsize Ah = 1/24. We start with the unstable uniform solution {u*,v*) with 
small random perturbation superimposed. Thus the initial profiles of u and 
V are completely random without any spatial correlation. And we perform a 
series of of two-dimensional simulations (figures [H [HlandUD, in each, the ini- 
tial condition was always a small amplitude random perturbation (±5 x 10~^) 
around the steady state {u*,v*), and patterns developed spontaneously. 

In the numerical simulations, different types of dynamics are observed and we 
have found that the distributions of predator and prey are always of the same 
type. Consequently, we can restrict our analysis of pattern formation to one 
distribution. In this section, we show the distribution of prey, for instance. 

From the analysis in section 2 and the bifurcation diagram (figure [T1(A)), 
the results of numerical simulations show that when parameters a, (3, c?2 are 
determined, the type of the system dynamics is determined by the values of 
7 and di. And for different sets of parameters, the features of the spatial 
patterns become essentially different when 7 exceeds the critical value 'Jh and 
7t respectively, which depend on di. 

Figure [21 shows the evolution of the spatial pattern of prey at 0, 10000, 20000 
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(C) (D) 

Fig. 2. Grey-scaled snapshots of contour pictures of the time evolution of the prey 
of model (g]) with 7 = 6.0 > 7//. (A) iteration; (B) 10000 iteration; (C) 20000 
iteration; (D) 50000 iteration. 

and 50000 iterations with 7 = 6.0, more than the Hopf bifurcation thresh- 
old ''^H = 5.552147102 and less than the Turing bifurcation threshold 7^ = 
7.265163898. In this case, one can see that for model (jl]), the random initial 
distribution around the steady state (m*, v*) = (0.39965, 0.26392) leads to the 
formation of the spiral wave pattern in the domain (figure [2](D)). In other 
words, in this situation, spatially uniform steady-state predator-prey coexis- 
tence is no longer. Small random fluctuations will be strongly amplified by 
diffusion, leading to nonuniform population distributions. From the analysis 
in section 2, we find with these parameters in domain II, the spiral pattern 
arises from Hopf instability. 

When 7 = 8.0 > 7t = 7.265163898, in this case, parameters in domain 
III (figure [U^A)), both Hopf and Turing instabilities occur. The nontrivial 
stationary state is {u*, v*) = (0.29974, 0.23088). As an example, the formation 
of a regular macroscopic two-dimensional spatial pattern, the chaotic spiral 
pattern, is shown in figure [31 
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(C) (D) 

Fig. 3. Grey-scaled snapshots of contour pictures of the time evohition of the prey 
of model 1^ with 7 = 8.0 > -yr- (A) iteration; (B) 10000 iteration; (C) 20000 
iteration; (D) 50000 iteration. 

Comparing this situation (figure [3]) with the one above (figure [2]), we can see 
that the formation of chaotic spiral patterns (figure [3]^C, D)) are caused by 
Turing instability. 

For the sake of learning the dynamics of model (jl]) further, we illustrate the 
phase portraits and time-series plots in figure |H From figure |1](A), for 7 = 6.0, 
one can see that a quasi limit cycle arises, which is caused by the Hopf bi- 
furcation. Furthermore, we can calculate that the frequency of periodic os- 
cillations in time (figure 111(B)) is u; = 27.5899, and corresponding wavelength 
A = 0.2277. And from ([8]) and iQ, we know that at the critical value of Hopf bi- 
furcation 'jH, the frequency of the periodic oscillations in time is uh = 24.0748, 
the corresponding wavelength Xh = 0.2610. When 7 = 8.0, the dynamical be- 
havior are shown in figures |1](C) and (D). From figure |1](C), one can see that 
there exhibits a "local" phase plane of the system invaded by the irregular 
spatiotemporal oscillations. Instead of the limit cycle in the case above (fig- 
ure |1](A)), as happens in the case of smooth pattern formation, the trajectory 
now fills nearly the whole domain inside the limit cycle. This regime of the 
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system dynamics corresponds to spatiotemporal chaos. And the spatial sym- 
metry of model (jl]) is broken and the patterns are oscillatory in space with 
the wavelength A = 0.1629 while at the critical value of Turing bifurcation, 
At = 0.2834. 




(C) (D) 

Fig. 4. Dynamical behavior of model ([4]). Phase planes of the model at a fixed point 
inside the domain occupied by irregular spatiotemporal oscillations with 7 = 6 (A) 
and 7 = 8 (C), respectively. Time-series plots with 7 = 6 (B) and 7 = 8 (D), 
respectively. 



The above results are obtained from the initial conditions which was always a 
small ar nphtude random perturbation around the steady state (u*, v*). In ref - 
erences ( iMedvinsky et all . l2002l : ISherratt. Eagan fc Lewisl . 119971 : iGarvid . 120071 ). 
the authors have studied the pattern formation arising from special initial 
conditions. They indicated that the spatiotemporal dynamics of a diffusion- 
reaction system depends on the choice of initial conditions. And the initial 
conditions are deliberately chosen to be asymme tric in order to rnake a ny in- 
fluence of the corners of the domain more visible (IMedvinsky et all . |2002| ). The 
initial localized introduction of predators into a uniform dist ribution of prey 
led to the spread of predators over the domain (IGarvid . 120071 ). An important 
new feature in the two-dimensional solutions is the way in which asymmetries 
in the initial introduction of predators are reflected in the long-time solutions. 



Based on the discussion above, we employ three categories of initial perturba- 
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tions (figure [5]) for furtlier learning tlie evolution of the spatial pattern of prey 
of model (jl]). In both parts of figure [5l predators were initially introduced in a 
spatially asymmetric manner. In the following, model (jl]) was solved numeri- 
cally with a time stepsize of At = 1 and space stepsize Ah = 1/3. 



sa 




3^ 




(A) 



Fig. 5. Three categories of initial perturbations corresponding with figures 6 — 9. 

In the first, predators are introduced in a small, localized region of the circle 
domain (Figure [5t^ A)): 



u{x,y,0) = 1.0, 



v{x,y,0) 



0.2 (x - 200)2 + (y - 200)2 < 400, 



otherwise 



(12) 



The second category of initial perturbation that we have used is the introduc- 
tion of predators along a line, which is otherwise in the steady state {u*,v*) 
(figure EKB)): 



u{x, y, 0) = u* + 0.005 ■ exp(^-(x - 200)^ - {y - 200) 
v{x, y, 0) = V* — 0.005 ■ exp \x — y\^ . 



(13) 



Thirdly, we employ the so-called pitchfork initial conditions (figure [5](C)): 



u{x, y, 0)=u* + 0.005 ■ exp (^-(x - 200)^ - (y - 200)^ j , 
v{x,y,0) = v* - 0.005 ■ exp(^-^\{x -y){x + y- 400) 1^ 



(14) 



The numerical simulations results of pattern formation of model (jlj) with 
above three categories of initial perturbations are shown in figures [6] (7 = 6.0) 
and [3 (7 = 8.0). 
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Fig. 6. Grey-scaled snapshots of contour pictures of the time evolution of the prey 
of system at different instants with < 7 = 6.0 < 7t. (A)(B)(C) are the three 
category of initial conditions corresponding to Figure (jH) and the iterations are: (i) 
3000; (ii) 5000; (iii) 10000; (iv) 50000. 

In figure [6]( A) , in the case of the first category of initial perturbation which 
predators are introduced in a small and localized region of the circle domain, 
one can see that after a symmetrical target pattern (figure E](A)(i)), it grows 
slightly and the spiral pattern (exterior) with target pattern (interior) emerges 
(figure [6t^A)(ii, iii)), finally with the appearance of spiral pattern in the whole 
domain (figure [6]( A) (iv)). And in figure [7|^A), 7 = 8.0, with the same initial 
conditions, a target pattern (figure [7]( A) (i)) emerges, the destruction of the 
target begins from the center, and leads to the formation of the spiral pattern 
(interior) with target pattern (exterior) (figure [7|( A) (ii, iii)), finally, the chaotic 
spatial pattern prevails the whole domain (figure [T]^ A) (iv)). Comparing Fig- 
ure [6](A) (iv) with figure [7]( A) (iv) , we find that Hopf instability leads to the 
formation of spiral patterns and the Turing instability destroys the spiral pat- 
tern and leads to the formation of chaotic spatial patterns. Moreover, in these 
two cases, the initial nonuniformity spreading outwards through the domain 
from the center provides additional evidence for spatiotemporal pattern. 

In figure [6]^B), with initial condition (fT3|) . bi-target pattern (figure [6](B)(i)) 
emerges, and finally, the phase waves appear in the whole domain (figure [6]^B) 
(iv)). And in figure [Tl^B), more differently, after a target pattern(figure [Tl^B) 
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Fig. 7. Grey-scaled snapshots of contour pictures of the time evolution of the prey 
of system at different instants with 7 = 8.0 > 7t. (A)(B) (i) 1500; (ii) 3000; (iii) 
5000; (iv) 50000. (C) (i) 1000; (ii) 3000; (iii) 5000; (iv) 50000. 

(i)), a chaotic spatial pattern occurs (figure [7](B) (iii, iv)). 

In the case of the third category initial perturbation (IT^ . one can see that 
figure El^C) and[7](C) follow similar scenario to the previous case, respectively 
(figures E](B) and[71^B)). The differences are that there are four-target patterns 
(figure [6t^C)(i) and [6t^C)(ii)) while in the previous case bi-target patterns oc- 
cur. Comparing these two cases with the initial perturbations defined by equa- 
tions f|T3l) and f|T^ . we find that in the beginning of evolution of the spatial 
pattern of prey, the special initial conditions have an effect on the formation 
of spatial patterns, though the effect is less and less with the more and more 
iterations. 

In order to make it more clearer, we show space-time plots in figures[8](7 = 6.0) 
and [9]^7 = 8.0). The method of space-time plots is that let y he a constant 
(here, y = 200, the center line of each snapshots), from each pattern snapshots, 
choose the line y = 200, and pile these lines in-time-order. The space-time plots 
show the evolution process of the prey u throughout time t and space x. On 
the other hand, we have found that the distributions of predator and prey are 
always of the same type. In this section, we show the distribution of prey, for 
instance. 
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Fig. 8. Space-time plots of variable 7 = 6.0. Other parameters are the same as those 
in Figure[2] and [6j 

In figure [HI 7 = 6.0, other parameters equal figures [2] and El with the three 
categories of initial perturbations (figure [5]), space-time plots at different times 
are shown. The first row in figure Et^A-C) displays the time evolution of the 
prey with the iterations from to 5000, while the second row (figure [Ht^D-F)) 
displays the time evolution of the prey with the iterations from 45000 to 50000. 
The three columns correspond to the three categories of initial perturbations 
(figure [5]). Form figure [HI and figure El one can see that when 7 = 6.0, in this 
case, the Hopf instability occurs, the circular initial condition creates spiral 
wave (figure E1(A) and figure [8]^D)), and the other two initial conditions create 
phase wave (figure E](B, C) and figure [8]^E, F)). Furthermore, in this case, the 
second and third initial conditions have the same effect on pattern formation 
(figure EKE, F)). 



16 






Fig. 9. Space-time plots of variable 7 = 8.0. Other parameters are the same as those 
in Figure [3] and [71 

When 7 = 8.0, other parameters equal figures [3] and [71 space-time plots at 
different times with the three categories of initial perturbations (figure [5]) are 
shown in figure [9l In this case, both Hopf and Turing instability occur. Fig- 
ure [9](A-C), displays the time evolution of the prey with the iterations from 
to 5000, while the second row (figure [9](D-F)) displays the time evolution 
of the prey with the iterations from 45000 to 50000. The three columns corre- 
spond to the three categories of initial perturbations (figure [5]) . Form figure [9] 
and figure [TJ one can see that all three initial conditions have the same effect 
on pattern formation (figure [Tl^iv) and figure [9](D-F)). 



In our previous work (IWang. Liu fc Jinl . |2007| ). we indicated that for the ratio 



dependent predator-prey system with Michaelis-Menten-type functional re- 
sponse, the case of predator-dependent, the stationary patterns did not depen- 
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dent on the initial conditions. One can see that, in the case of prey-dependent 
type (jlj), pattern formation do depend on the initial conditions. This is a ma- 
jor difference between prey- dependent and predator-dependent predator-prey 
models. 



4 Discussions 



In this paper, we have presented a theoretical analysis of evolutionary pro- 
cesses that involve organisms distribution and their interaction of spatially 
distributed population with local diffusion. And the numerical simulations 
were consistent with the predictions drawn from the bifurcation analysis, that 
is, Hopf bifurcation and Turing bifurcation. 



If the parameter 7 in the domain II, the Hopf instability occurs, the destruction 
of the pattern begins from the interior, while it begins from the exterior if 7 
located in the domain I II, both Hopf and Turing instabilities occur. From 
an ecological viewpoint (ISherratt. Eagan fc Lewis! . Il997l ). it shows that the 
initial and relatively rapid invasion of prey by predators can be followed by 
two subsequent invasions. 

The effort to explain the distribution of populations in terms of the move- 
ments of individuals is an extension of one of the most successful applications 
of mathematics to ecological phenomena, the use of diffusion model to describe 
dispersal. The basic idea is that, although organisms do not move randomly, 
the collective behavior of large numbers of such individuals nmy_be indistin- 
guishable ( at the scale of the population) from what they did (Levin . 19921). In 
refere nces (lLevinl . ll992l : ISavill fc Hogewegi . ll999l : lAlonso. Bartumeus fc Catalan 
2OO2I ). the authors indicated that the basic idea of diffusion- driven instabil- 
ity in a reaction-diffusion system can be understood in terms of an activator 



inhib itor system. The functioning of this mechanism is based on three points (lAlonso. Bartumeus fc Cat 
2OO2I ). First, a random increase of activator species should have a positive ef- 
fect on the creation rate of both activator and inhibitor species. Second, an 
increment in inhibitor species should have a negative effect on formation rate 
of both species. Finally, inhibitor species must diffuse faster than activator 
species. Certainly, the reaction-diffusion predator-prey model (jl]), with Ivlev 
functional response and predators diffusing faster than prey, provides this 
mechanism. 



Spirals and curves are the most fascinating clusters to emerge from the predator- 
prey model. A spiral will form from a wave front when the rabbit line (which is 
leading the front) overlaps the pursuing line of predator. The prey on the ex- 
treme end of the line stop moving as there are no predator in their immediate 
vicinity. However the prey and the predator in the center of the line continue 
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moving forward. This forms a small trail of prey at one (or both) ends of the 
front. These prey start breeding and the trailing line of prey thickens and 
attracts the attention of predator at the end of the fox line that turn towards 
this new source of prey. Thus a spiral forms with predator on the inside and 
prey on the outside. If the original overlap of prey occurs at both ends of the 
line a double spiral will f orm. Spirals can also form as a p rey blob collapses 



after predator eat into it (iHawick. James fc Scogingd . l2006l ) . 



And a random increase of activator species (prey, u) has a positive effect on 
the creation rate of both activator and inhibitor species. Random fluctuations 
may cause a nonuniform prey density. This elevated prey density has a positive 
effect both on prey and predator population growth rates. From model ([3]), 
we can obtain per capita rates: 



if = (1 - t.) - ^(1 - e--"), if = /?(«- 1 - ae-^n- (15) 



Since the first equation in f[T^ is a one-humped function of prey density u, 
prey growth rate can be increased by a higher local prey density at least in 
a range of parameter values. The second equation in (fTSll . predator numerical 
response, is an ever-increasing function of u, and high prey density always has 
a positive influence on predator growth. More importantly, inhibitor species 
(predator, v) must diffuse faster than activator species (prey, u), for an in- 
crement in inhibitor species may have a negative effect on formation rate of 
both species. Thus, as random fluctuations increase local prey density over its 
equilibrium value, prey population undergoes an accelerated growth. Simul- 
taneously, predator population also increases, but as predators diffuse faster 
than prey, they disperse away from the center of prey outbreaks. If relative 
diffusion (^2/^1) is large enough, prey growth rate will reach negative values 
and prey population will be driven by predators to a very low level in those 
regions. In other words, where the prey density is at their maximal value dif- 
fusion will lower the prey density at that point. Conversely, where the prey 
density is at their minimal value diffusion will increase the prey density at that 
point. That is, prey flow from high density to low density regions in space. 
Moreover, the faster the diffusion the greater the flow. When the prey density 
is high (in fact, when V^m < 0) proportionately less of the slower diffusing 
prey leaves these points in space than the faster diffusion prey and, therefore, 
the proportion of the slower diffusing prey at that point increases more than 
the proportion of the faster diffusing prey. Conversely, at low prey density 
(when V^M > 0) proportionately more faster diffusing prey enter a point in 
space and the proportion of these prey increases more than the proportion of 
the slower diffusing prey. Hence, at a given position in space, when the prey 
density is high the proportion of slower diffusing prey increases and when the 
prey density is low the proportion of the faster diffusing prey increases. There- 
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fore we see oscillations in the proportions of the prey at the same frequency 
as the oscillations of the density waves. The final result is the formation of 
patches of high prey density is surrounded by areas of low prey density. And 
predators follow the same pattern. 



On the other hand, in two-dimensional reaction-diffusion systems, rotationally 
symmetric patterns, known as targets or sinks, and a generalization of them 
with broken circular symmetry, spirals are being investigated experimentally 
as well as theoretically in many nonlinear systems. The Belousov-Zabotinsky 
reaction is a well investigated excitable reaction-diffusion system that shows all 
these patterns. Spirals are characteristic patterns in slime mold aggregates and 
are an important observation in cardiac arrythmias as well. Targets and spi- 
rals, which are generally found to form around some defects, p recede some de- 
fect m ediated chaos, commonly known as spiral defect chaos (jBhattacharyayl . 
200lh . 



Particularly, spiral patterns are being investigated theoretically in a number of 

reaction-diffusion predator -prey system s, such as HoUing-type model (ISavill &: Hogewed. 
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1999 : lMalchow et all.l2000[). lylev-tvpe (jSherratt, Lewis &: Fowlerl.ll995l:ISherratt. Eagan fc Lewis 



Kav fc Sherrattl . jgOOd : IPearce et ai I2OO6I : ICarviel . l2007l : IPreedv et al 



20061 : lUriu fc Iwasal . l2007h . and so on. The functional responses of these predator- 
prey models are all prey-dependent. It is necessary or coincident? 

It is well known that, for react ion- diffusion predator-prey systems, under 
suitable conditions, the destabilized uniform distributions give way to stable 

nonuniform patterns, whic h can provid e the local inf ormation of that specifies 

patte rns of differentiation (jLevinl . ll992l ). In reference (lAlonso. Bartumeus fc Catalan! . 
20021 ). the authors also indicted that a simple general model for predator-prey 
dynamics with predator-dependent functional response, a reaction-diffusion 
system that could develop diffusion- driven instabilities. On the contrary, if 
the functional response depends only on prey density, diffusion instabilities 
are not possible. In fact, in the case of predator-dependent, the interactions 
between dispersing populations induce spatial heterogeneity and/or tempo- 
ral fluctuations through so-called self-structuring without help from external 
forcing, and the patterns are endogenous. While the prey- dependent predator- 
prey models are self-oscillation ones which are called oscillatory systems, there 
are typical patterns including spiral waves, turbulence, and target patterns. 
This is the reason why we can find the spiral and target waves in model (jl]). 



From the theore t ical st udy on the three-dimensional patterns (iLeppanen et all . 
2OO3I : iLeppanenl . |200J) , it is possible that the three-dimensional patterns can 
be reflected by the two-dimensional patterns. And two-dimensional patterns 
might be sufficient to understand the g e neral properties of dissipative struc- 
tures (IShoji. Yamada. Ueyama fc Ohtal . 120071 1. So our two-dimensional spa- 
tial patterns may indicate the vital role of pattern formation in the three- 
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dimensional spatiotemporal organization of the predator-prey system. 



Furth ermore, from references (jSherratt. Perumpanani fc Oweru . ll999l : lByrne et al 



20061 ). we think that pattern formation of spatial model (jlj) with special choice 
initial conditions (IT^ - (IT^ can be used to explain other diffusion process, such 
as tumor growth, and so on. 

On the other hand, the ecosystem is so complicated that we cannot use a 
single method to study. We must use mixed methods, such as analytical or 
experimental or numerical method. 
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